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A CONTOUR-INTEGRAL BASED QZ ALGORITHM FOR 
GENERALIZED EIGENVALUE PROBLEMS 

GUOJIAN YIN* 


Abstract. Recently, a kind of eigensolvers based on contour integral were developed for com¬ 
puting the eigenvalues inside a given region in the complex plane. The CIRR method is a classic 
example among this kind of methods. In this paper, we propose a contour-integral based QZ method 
which is also devoted to computing partial spectrum of generalized eigenvalue problems. Our new 
method takes advantage of the technique in the CIRR method of constructing a particular subspace 
containing the eigenspace of interest via contour integrals. The main difference between our method 
and CIRR is the mechanism of extracting the desired eigenpairs. We establish the related framework 
and address some implementation issues so as to make the resulting method applicable in practical 
implementations. Numerical experiments are reported to illustrate the numerical performance of our 
new method. 
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1. Introduction. Let A and B be large n x n matrices. Assume that we have 
a generalized eigenvalue problem 

Ax = XBx, (1-1) 

and want to compute the eigenvalues Xi, along with their eigenvectors x^, of (HU) 
inside a given region in the complex plane. This problem arises in various areas of 
scientihc and engineering applications, for example in the model reduction of a linear 
dynamical system, one needs to know the response over a range of frequencies, see 
[SlIISlET]. Computing a number of interior eigenvalues of a large problem remains one 
of the most difficult problems in computational linear algebra today [10] . In practice, 
the methods of choice are always based on the projection techniques, the key to the 
success of which is to construct an approximately invariant subspace enclosing the 
eigenspace of interest. The Krylov subspace methods in conjunction with the spectral 
transformation techniques, such as the shift-and-invert technique, are most often used 
[221 [26]. 

Recently, the eigensolvers based on contour integral were developed to compute 
the eigenvalues inside a prescribed domain in the complex plane. The best-known 
methods of this kind are the Sakurai-Sugiura (SS) method [24] and the FEAST algo¬ 
rithm [20] . A major computational advantage of these contour-integral based meth¬ 
ods is that they can be easily implemented in modern distributed parallel computers 
mm- The FEAST algorithm works under the conditions that matrices A and B are 
Hermitian and B is positive definite. In the SS method, the original eigenproblem 
(HU is reduced to a small one with Hankel matrices, if the number of sought-after 
eigenvalues is small. However, since Hankel matrices are usually ill-conditioned [5], 
the SS method always suffers from numerical instability mm- By noticing this fact, 
later in [25], Sakurai et al. used the Rayleigh-Ritz procedure to replace the Hankel 
matrix approach to get a more stable algorithm called CIRR. 
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Originally, the CIRR method was formulated under the assumptions that matrices 
A and R are Hermitian and B is positive dehnite, i.e., (11.11) is a Hermitian problem 
[3]. Moreover, it is required that the eigenvalues of interest are distinct. In [T7], the 
authors adapted the CIRR method to non-Hermitian cases; meanwhile, they presented 
a block version of the CIRR method so as to deal with the degenerate systems. 

The CIRR method is always accurate and powerful. It first constructs a sub¬ 
space containing the eigenspace of interest through a sequence of particular contour 
integrals. Then the orthogonal projection technique is used to extract desired eigen- 
pairs. In our work, we propose a contour-integral based QZ method for solving partial 
spectrum of dni). The motivation stems from the attempt of using the oblique pro¬ 
jection method, instead of the orthogonal one, to extract desired eigenpairs in the 
CIRR method. When using oblique projection technique, the most important task is 
to find an appropriate left subspace, we borrow ideas of the JDQZ method m, and 
derive our new method. We establish the related mathematical framework. Some 
implementation issues will also be discussed before giving the resulting algorithm. 

The rest of the paper is organized as follows. In Section 2, we briefly review 
the CIRR method [25]. In Section 3, we derive a contour-integral based QZ method 
and establish the related mathematical framework. Then we will discuss some im¬ 
plementation issues and present the complete algorithm. Numerical experiments are 
reported in Section 4 to illustrate the numerical performance of our new method. 

Throughout the paper, we use the following notation and terminology. The sub¬ 
space spanned by the columns of matrix X is denoted by span{X}. The rank of 
matrix A is denoted by rank(A). For any matrix S, we denote the submatrix that lies 
in the first i rows and the first j columns of S by the submatrix consisting 

of the first j columns of S by and the submatrix consisting of the first i rows 

of S by The algorithms are presented in Matlab style. 

2. The CIRR method. In [24], Sakurai et al. used a moment-based technique 
to formulate a contour-integral based method, i.e., the SS method, for finding the 
eigenvalues of cu inside a given region. In order to improve the numerical stability 
of the SS method, a variant of it used the Rayleigh-Ritz procedure to extract desired 
eigenpairs. This leads to the so-called CIRR method [TTI [25] . Originally the CIRR 
method was derived in |25j under the assumptions that (i) matrices A and B are Her¬ 
mitian with B being positive definite, and (ii) the eigenvalues inside the given region 
are distinct. In the authors adapted the CIRR method to the non-Hermitian 
cases, meanwhile, a block version was proposed to deal with the degenerate problems. 
In this section we give a briefly review of the block CIRR method. 

The matrix pencil zB — H is regular if det( 2 :i? — A) is not identically zero for 
all z £ C [2l|9] . The Weierstrass canonical form of regular matrix pencil zB — H is 
defined as follows. 

Theorem 2.1 m )- Let zB — A be a regular matrix pencil of order n. Then 
there exist nonsingular matrices S and T £ such that 


TAS = 


Jd 

0 


0 

h^n—d 


and 


TBS = 


Id 

0 


0 

Nn-d 


( 2 . 1 ) 


where Jd is a d x d matrix in Jordan canonical form with its diagonal entries corre¬ 
sponding to the eigenvalues of zB — A, Nn-d is an (n — d) x (n — d) nilpotent matrix 
also in Jordan canonical form, and Id denotes the identity matrix of order d. 
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Let Jd in (EU be of the form 


Jd 


Jdi(Ai) 0 

0 Jd2i^2) 


0 

0 


0 0 ••• Jdmi^m) 


( 2 . 2 ) 


where di = d, 1 < di < d for i = 1,... m and Jd^ (Ai) are di x di matrices of the 
form 


JdA^r) 


1 0 

0 A, 1 


0 ■ 

0 ’ 


: 1 
0 • •• 0 A* 


z = 1, 2,..., m 


with Xi being the eigenvalues. Here the Xi are not necessarily distinct and can be 
repeated according to their multiplicities. 

Let us partition S into block form 


S=[S^,S2, 


Sm,S, 


m+lj ^ 


(2.3) 


where Si £ 1 < i < m, and S'm+i £ Then the first column in each 

Si is an eigenvector associated with eigenvalue Xi for z = 1,..., m [11I711I1117]. 

Let r be a given positively oriented simple closed curve in the complex plane. 
Below we show how to use the block CIRR method to compute the eigenvalues of 
(El) inside L, along with their associated eigenvectors. Without loss of generality, let 
the set of eigenvalues of El enclosed by L be {Ai,..., A;}, and s := di + ^2 H— ■ + di 
be the number of eigenvalues inside L with multiplicity taken into account. 

Define the contour integrals 


Fk := 


^{zB — A) ^Bdz, fc = 0,1,_ 


(2.4) 


2Tr^/^jr 

With the help of residue theorem in complex analysis [1], it was shown in m that 
= 5(,i,«)(J(i,«,i,,))'=(5-i)(i,,,), fc = 0,l,.... (2.5) 


Let h and g be two positive integers satisfying hg ^ s, and T be an n x h random 
matrix. Define 


Uk:=FkY, k = 0,...,g-l, and C/:= [C/q, Di,..., C/g_i]. (2.6) 

We have the following result for the CIRR method. 

Theorem 2.2. Let the eigenvalues inside T be Xi,Xi, then the number of 
eigenvalues of EP inside T is s, counting multiplicity. //rank([/) = s, then we have 

span{C/} = span{S'(._i.s)}. (2.7) 
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Proof. By (12.51) and (12.61) . we know that 

U = (2.8) 

where 

E = [iS-^)^i,s,)Y, Jil..s,l:s)iS-^)il:s,)Y, . . . , () ( 1 :.,) F] . (2.9) 

Since the rank of C/ is s, we have that E is full-rank, following from which the expres¬ 
sion (12.71) holds. □ 

According to Theorem 12.21 we know that spanjU} contains the eigenspace cor¬ 
responding to the desired eigenvalues. The block CIRR method uses the well-known 
orthogonal projection technique to extract the eigenpairs inside T from spanjU}, i.e., 
imposing the Ritz-Galerkin condition: 

Ax — XBx T span{17}, (2-10) 


where A € C and x € spanjU}. 

The main task of the block CIRR method is to evaluate Uk (cf. (12.61) 1. In practice, 
Uk have to be computed approximately by a numerical integration scheme: 

I 9 

= fc = 0,l,...,g-I, (2.11) 

^ 1=1 

where Zj are the integration points and ojj are the corresponding weights. From (I2.11|) . 
it is easy to see that the dominant work of the block CIRR method is actually solving 
q generalized shifted linear systems of the form 

{ZjB - A)X, = BY, j = l,2,...,g. (2.12) 

Noticing that the integration points Zj and the columns of right-hand sides are inde¬ 
pendent, the CIRR method can be easily implemented in modern distributed parallel 
computer. 

The complete block CIRR method is summarized as follows. 

Algorithm 1: The block CIRR method 
Input: h,g,q,Y G 

Output: Approximate eigenpairs (Ai,Xi), Ai inside T. 

1. Compute Uk,k = 0,1,g — 1, approximately by (12.111) . 

2. Compute the singular value decomposition of 17 = [Co,..., Cg-i] : U = UYV. 

3. Set A = U*AU and B = U*BU. 

4. Solve the generalized eigenproblem of size hg: Ay = \By, to obtain the 
eigenpairs {(Ai,yi)}Jfi. 

5. Compute Xi = Uyi, and select s approximate eigenpairs (Ai,Xi). 

3. A contour-integral based QZ algorithm. The contour-integral based 
methods are recent efforts for the eigenvalue problems. The CIRR method is a typical 
example among the methods of this kind. According to the brief description in the 
previous section, the basic idea of the block CIRR method can be summarized as 
follows: (i) constructing a particular subspace that contains the desired eigenspace 
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by means of a sequence of contour integrals (cf. (12.41) 1. and (ii) using the orthogonal 
projection technique, with respect to the subspace spanjU} (cf. (12.61) 1. to extract 
the desired eigenpairs. In this section, we will derive another contour-integral based 
eigensolver. The idea stems from the attempt to use the oblique projection technique 
to extract desired eigenvalues in the block CIRR method. Applying the oblique pro¬ 
jection method, the key step is finding a suitable left subspace. We find an appropriate 
left subspace via using the QZ method to generate a generalized Schur decomposition 
associated with the desired eigenvalues. This intention finally leads us to a contour- 
integral based QZ method for solving dm. We call the resulting algorithm CIQZ for 
ease of reference. 

In this section, we first detail the derivation of our contour-integral based QZ 
method. Later on, we discuss some implementation issues that our contour-integral 
based QZ method may encounter in the practical application, after that, we give the 
complete CIQZ method. 

3.1. The derivation of the CIQZ algorithm. The CIRR method uses the 
orthogonal projection technique to extract the sought-after eigenpairs from span{t/}. 
Here we consider using the oblique projection technique [4l[22], another class of pro¬ 
jection method, to compute the desired eigenpairs. 

Since spanjC/} contains the eigenspace of interest, it is natural to choose spanjC} 
as the right subspace (or search subspace). The oblique projection technique extracts 
the desired eigenpairs from span{{7} by imposing the Petrov-Galerkin condition, 
which requires orthogonality with respect to some left subspace (or test subspace)., 
say, spanjlT}: 


Ax — ARx T span{IT}, (3-1) 

where A is located inside T, x G spanjC/j, and IT is an n x s orthogonal matrix. 
Let T be an n X s matrix whose columns form an orthogonal basis of span{f/}. The 
orthogonality condition dm leads to the projected eigenproblem 

W*AVy = XW*BVy, (3.2) 


where y £ C® satisfies x = Ty. 

Now our task is to seek an appropriate left subspace span{IT}. Our discussion 
begins with a partial generalized Schur form for matrix pair {A, B). 

Definition 3.1 m )- A partial generalized Schur form of dimension s for a 
matrix pair {A, B) is the decomposition 

AQs = Z,Hs, BQs = Z,Gs, (3.3) 

where Qs and Zg are orthogonal n x s matrices, and Hg and Gg are upper triangular 
s X s matrices. A column {Qs){-.,i) is referred to as a generalized Schur vector, and we 
refer to a pair ((Qs)(:,i), (lls)(i,i)/(Gs)(i^i)) as a generalized Schur pair. 

The formulation (13.31) is equivalent to 

{Zg)*AQg=Hg, {Zg)*BQg=Gg, (3.4) 

from which we know that {Hs)(^i^i)/{Gs)(^i^i) are the eigenvalues of {Hs,Gg). Let yi 
be the eigenvectors of pair {Hg,Gg) associated with (iL;,)(j i)/(Gs)(i_i), then we have 
((-ffs)(i,i)/(Gs)(i,q,(5syi) are the eigenpairs of {A,B) [TTIIIS- 
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Applying the QZ algorithm to (13.21) to yield generalized Schur form 

{Pl)*{W*AV)Pr = Ha and {Pr)* iW* BV)Pr = Hr, (3.5) 

where Pr and Pr are orthogonal s x s matrices, Ha and Hr are upper triangular 
s X s matrices. The eigenvalues of pair {W*AV,W*BV) are 

[lais]. 

Comparing (13.4|) with (13.51) . it is readily to see that we have constructed a partial 
generalized Schur form in (13.51) for matrix pair {A,B): VPr constructs a Qg and 
WPr constructs a Zg. 

Since the desired eigenvalues are finite, the diagonal entries of Ha and Hr are 
non-zero, which means that Ha and Hr are nonsingular. In view of (j3.5l) . we can 
conclude that 


span{ITPL} = spanjACPR} = spa.n{BVPR}. (3.6) 

On the other hand, since Pr and Pr are nonsingular, we have 

span{IT} = spanjAId} = span{i3C}. (3.7) 

Motivated by (1X71) . we choose the left subspace spanjlT} to be spanjAU -f BU}. 
Below we want to justify this choice. 

Theorem 3.2. Let L,D G > s, be arbitrary matrices, and R = FqD. A 

projected matrix pencil zB — A is defined by B = L*BR and A = L*AR. If ranks 
of both T*(T“^)(,_i,s) and (>S'“^)(i:s^:)I? are s, then the eigenvalues of zB — A are 
Xi,..., Xi, i.e., the eigenvalues that are located inside T. 

The proof is almost identical with that of Theorem 4 in where the contour 
integrals Fk were defined as j>^ z^(zB—A)~^dz, that is, the term B was dropped 

comparing with the expression (12.51) . 

Theorem l3.2l savs that the desired eigenvalues can be solved via computing 

the eigenvalues of projected eigenproblem zB — A, if the ranks of both T*(T“^)(,_i,s) 
and (>S'“^)(i:s,:)Tl are s. Due to this, we want to show the following results. 

Theorem 3.3. If the rank of U is s, then the ranks of {AU + BU)*{T~^)(^,A:g) 
and (>S'“^)(i:s_:)C/ are s. 

Proof We first show that the rank of (*S'“^)(i:s_:)C/ is s. By (12.11) and (12.8L we 
have 

{S-^)il:s,)U = (5-')(l:s,:)5(,l:s)i5 = E. (3.8) 

Since U is full-rank, by (12.81) . we know that rank(i?) = s. Therefore, the rank of 

('S'"^)(l:s.:)C7 is S. 

Next we show that the rank of {AU -I- B{7)*(T“^)(. is s. For convenience, we 
turn to show that the rank of {{T~^)(^,A-.g))*{AU + BU), i.e., the conjugate transpose 
oi {AU + BU)*{T-^)(^,a:s)As s. 

Since spanjAU} = span{i3[/} (cf. (13.71) 1. there exists a hg x hg nonsingular 
matrix A such that AU = BU A. According to (12.IL (12.5|) . and (12.81) . we have 

{{T-^)i:,i,g))*{AU + BU) = {BS(^,A-.s)yBS^..A..g)E{A + R). (3.9) 

In view of (j2.1L we know BS(^,a.s) is full-rank, which means {BS(^,a.s))* ES(^:A:s'j is 
nonsingular. By (13.9p . we can conclude that (T“^)(. i.s))*(A[/ -I- BU) is full rank, 
thus the rank of {AU + B[/)*(T“^)(. is s. □ 
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Based on Theorem 13.21 and Theorem l3.31 we have that the eigenvalues of {{AU + 
BU)*AU, {AU + BU)*BU) are the eigenvalues of (11.11) inside F, which justifies our 
choice of taking the left subspace to be span{At/ -\- BU}. On the other hand, the 
columns of V and W form the base of span{AC/ + BU} and span{t7}, respectively. 
As a consequence, there exist hg x hg nonsingular matrices Pi and P 2 , such that 

[AU + BU)Pi = [IT, 0], UP 2 = [F, 0]. (3.10) 


Now, we have 


P*{z{AU + BU)*BU - {AU + BU)*AU)P 2 


zW*BV - W*AV O' 

0 0 


(3.11) 


Therefore, {W*AV, W*BV) shares the same eigenvalues with {{AU + BU)*AU, {AU + 
BU)*BU), which are {{HA){i,z)/{HB)(i,i)}t=i by dSH). Let {{HA){i,i)/{HB){i,i),yi) 
be the eigenpairs of {Ha,Hb), then according to (13.411 and (13.5L we have that 
{{HA){i,i)/{HB)(i,i),VPRyi) are exactly the eigenpairs of (II. 1|) inside F. 

We use the following algorithm to summarize the above discussion. 

Algorithm 2: A contour-integral based QZ algorithm. 

Input: h,g,q,Y G 

Output: Approximate eigenpairs (A^, x^), i = 1,..., s. 

1. Compute Uk, k = 0,1,..., g — 1, approximately by (12.111) . 

2. Form U = \U{j, Ui ,..., C/g-i] and compute orthogonalization: 

V = orth(t7) and W = orth(AlA -|- BV). 

3. Compute A = W*AV and B = W*BV. 

4. Compute [Sa, Sb,Ul, Ur, Vl,Vr] = qz(A, B). 

5. Compute Xi = {SA){i,i)/{SB){i,i) and = VUR{VR)(^.,^i). 

6. Select the approximate eigenpairs (Ai,Xi). 

3.2. The implementation issues. If we apply Algorithm 2 to compute the 
eigenvalues inside F, we will encounter some issues in practical implementation, just 
like other contour-integral based eigensolvers [20l[2l[25]. In this section, we discuss 
the implementation issues of our new method. 

The first issue we have to treat is about selecting a suitable size for the starting 
matrix Y, with a prescribed parameter g. Since U (cf. 12.6|) is expected to span a 
subspace that contains the eigenspace of interest, we have to choose a parameter h, 
the number of columns of Y, such that hg ^ s, the number of eigenvalues inside F. A 
strategy was proposed in [23] for finding a suitable parameter h for the block CIRR 
method. It starts with finding an estimation to s. Giving a positive integer /ip, by 
“Y/jp ~ N(0,1)”, we mean is an n x /ip matrix with i.i.d. entries drawn from 
standard normal distribution N(0,1). By (12.51) and (12.6L one can easily verify that 
the mean 


E[trace((F/,p)*FpY)ip)] = hp • trace(Fp) = hp • trace(S'(,_i,s)(S' ^)(i:s,:)) = /ip-s. (3.12) 
Therefore, 


sp := — • E[trace((Yho)*FpYjip)] 

riQ 
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(3.13) 









gives an initial estimation to s [I1I27]. With this information on hand, the strategy 
in works as follows: (i) set h = 1"-^], where k > 1, (ii) select the starting matrix 
Y G and compute Uk by (12.lip , (hi) if the minimum singular value amin of 

U = [C/q: ■ • ■ jUg-i] is small enough, we find a suitable h; otherwise, replace h with 
ah and repeat (ii) and (iii). We observe that the formula (13.131) always gives a good 
estimation of s. However the computed Sq may be much larger than s in some cases, 
such as the matrices A and B are ill-conditioned, which leads to that it is potentially 
expensive to compute the singular value decomposition of U. Due to this fact, in our 
method we turn to use the strategy proposed in m, whose working mechanism is as 
follows: use the rank-revealing QR factorization mm to monitor the numerical rank 
of {7, if U is numerically rank-deficient, then it means that the subspace spanned by 
U already contains the desired eigenspace sufficiently, as a result, we find a suitable 
parameter h. 

Another issue we have to address is designing the stopping criteria. The stopping 
criteria here include two aspects: (i) all computed approximate eigenpairs attain the 
prescribed accuracy, and (ii) all eigenpairs inside the given region are found. 

As for the first aspect of the stopping criteria, since we can only compute U 
approximately by some quadrature scheme (cf. (I2.11[) l. the approximate eigenpairs 
computed by Algorithm 2 may be unable to attain the prescribed accuracy in practical 
applications. A natural solution is to refine U (step 2 in Algorithm 2) iteratively. A 
refinement scheme was suggested in [16]. Let = Y and I be a positive integer, 
the refinement scheme iteratively computes Uj:'^ = FkUQ by a q-point numerical 
integration scheme: 

= fc = 0,l,...,5-l, (3.14) 


and then constructs 


= 


7j(0 jjil) 


,u. 


(0 


s-1 


(3.15) 


The refined is used to form projected eigenproblem (13.2p . through which we 
compute the approximate eigenpairs. The accuracy of approximate eigenpairs will be 
improved as the iterations proceed, see [23] for more details. 

If all s approximate eigenpairs attain the prescribed accuracy after a certain 
iteration, we could stop the iteration process. However, in general we do not know 
the number of eigenvalues inside the target region in advance. This fact leads to the 
second aspect of the stopping criteria: how to guarantee that all desired eigenpairs are 
found when the iteration process stops. We take advantage of the idea proposed in 
m- The rationale of the idea is that, as the iteration process proceeds, the accuracy 
of desired eigenpairs will be improved while the spurious ones do not, as a result, 
there will exist a gap of accuracy between the desired eigenpairs and the spurious 
ones m- Based on this observation, a test tolerance ry, say 1.0 x 10 is introduced 
to discriminate between the desired eigenpairs and the spurious ones. Specifically, for 
approximate eigenpair (Ai,Xi), define the corresponding residual norm as 


IIAx,; - XjBXiW 

IIAxill -f IlHxill’ 


(3.16) 


If ri < ry, then we view (Ai,Xi) as an approximation to a sought-after eigenpair and 
refer to it as a filtered eigenpair by rj. If the numbers of filtered eigenpairs are the 












same in two consecutive iterations, then we set them to be the number of eigenvalues 
inside F, see m for more details. 

From (I3.14|) we can see that, in each iteration, the dominate work is to compute 
q generalized shifted linear systems of the form 

iz,B - ^=l,2,...,g. (3.17) 


Integrating the above strategies with Algorithm 2, below we give the complete 
CIQZ algorithm for computing the eigenpairs inside the given region F. 

Algorithm 3: The complete CIQZ method 

Input: A, B, ho, g, q, K,r],e, maxAtev. 

Output: Approximate eigenpairs (Aj, x^), i = 1,..., s. 

1. Let Yhg ^ N(0,1), compute t/fe, fc = 0,..., g — 1, by (12.111) . 

2. Compute sq = |"j^trace((y/j(,)*i/^o)], and set h = max{, ho}. 

3. If h > ho 

4. Pick Yh-ho ~ N(0,1) and compute Uk by by (12.lip . Augment Uk 


Uk,Uk 


£ ami construct U = 

Uq, Cl, ... , Ug-l 


Uq, Cl, ... , Ug-l 


to Uk- Uk = 

5. Else 

6. Set h = ho and construct C = 

7. End 

8. Compute the rank-revealing QR factorization: C = VRIV. Set si = rank(i?). 
If Si < hg, stop; otherwise, set ho = h, h = nh and go to step 3. 

9. Set e(0) = 0 and c(0) = n. 

10. For fc = 1, 2,..., maxJter 

11. Compute the orthogonalization: V = orth(C) and W = orth(AE -|- BV). 

12. Compute A = W*AV and B = W*BV. Set si = rank(A). 

13. Compute [Sa, Sb, Ul, Ur, Vl,Vr] = qz(A, B). 

14. Compute A = and A = ECij(VR)(._i),i = l,...,Si. 

15. Set r = [ ], A^^^ = [ ], and c(fc) = 0. 

16. For z = 1 : Si 

17. Compute u = ||Ax,: - ACxi||/(||Ax,:|| + ||Sxi||)- 

18. If Xi inside F and < g, then c{k) = c{k) -|- 1, r = [r, ri], 
x(k) ^ [xW,Xi] and AW = [A^, A]. 

19. End 

20. Set e{k) = max(r). 

21. If c{k) = c{k— 1) and e{k) < e, set Xi = Stop. 

22. Set Y = Uo, and compute Uk by (12.111) . Construct C = Co, Ui,..., Cg-i 

23. End 


Here we give some remarks on Algorithm 3. 

1. Steps 1 to 8 are devoted to determining a suitable parameter h for the starting 
matrix Y. Meanwhile, a matrix C is also generated. 

2. The for-loop, steps 16 to 19, is used to detect the spurious eigenvalues. Only 
the approximate eigenpairs whose residual norms are less than g are retained. 

3. Step 21 refers to the stopping criteria, which contain two aspects: (i) the 
number of filtered eigenpairs by g is the same with the one in the previous 
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Table 4.1 

Test problems from Matrix Market that are used in our experiments. 


No. 

Problem Size 

Matrix nnz Property cond 

1 

BFW398 398 

A 3678 real unsymmetric 7.58 x 

B 2910 real symmetric indefinite 3.64 x 10^ 

2 

BFW782 782 

A 7514 real unsymmetric 4.63 x 10^ 

B 5982 real symmetric indefinite 3.05 x 10^ 

3 

DWG961 961 

A 3405 complex symmetric indefinite Inf 

B 10591 complex symmetric indefinite 3.21 x 10^ 

4 

MHD1280 1280 

A 47906 complex unsymmetric 9.97 X 10^^ 

B 22778 complex Hermitian 5.99 X 10^^ 

5 

MHD3200 3200 

A 68026 real unsymmetric 2.02 X 10‘‘‘‘ 

B 18316 real symmetric indefinite 2.02 X 10^® 

6 

MHD4800 4800 

A 102252 real unsymmetric 2.54 X lO'^' 

B 27520 real symmetric indefinite 1.03 X 10^^ 


iteration, and (ii) the residual norms of all filtered eigenpairs are less than 
the prescribed tolerance e. 

4. By (13.7L theoretically, the left subspace can be chosen either span{A?7} or 
span{i?17}. However, in practical implementation U can only be computed 
by a quadrature scheme to get an approximation U, in Algorithm 3 we choose 
the left subspace to be span{Af/ + Ht/} so as to include the information of 
both span{At/} and span{i?H}. 

4. Numerical Experiments. In this section, we use some numerical exper¬ 
iments to illustrate the performance of our CIQZ method (Algorithm 3). The test 
problems are from the Matrix Market collection . They are the real-world problems 
from scientific and engineering applications. The descriptions of the related matrices 
are presented in Table 14.11 where nnz denotes the number of non-zero entries and 
cond denotes the condition numbers which are computed by Matlab function condest. 
All computations are carried out in Matlab version R2014b on a MacBook with an 
Intel Core i5 2.5 GHz processor and 8 GB RAM. 

We use Gauss-Legendre quadrature rule with q = 16 quadrature points on T 
to compute the contour integrals (I3.14|) [ 8 ]. As for solving the generalized shifted 
linear systems of the form (I3.17|) . we first use the Matlab function lu to compute 
the LU decomposition of A — ZjB,j = 1, 2,..., g, and then perform the triangular 
substitutions to get the corresponding solutions. In the experiments, the size of 
sampling vectors /iq a-nd the parameter g are taken to be 20 and 5, respectively. 

Experiment 4.1 The goal of this experiment is to show the convergence behavior 
of CIQZ. The test problem is the bounded fineline dielectric waveguide generalized 
eigenproblem BFW782 (cf. Table |4T]) [6]. It stems from a finite element discretization 
of the Maxwell equation for propagating modes and magnetic field profiles of a rectan¬ 
gular waveguide filled with dielectric and PEC structures [4]. We are interested in the 
eigenvalues inside the circle T with center at 7 = — 6.0 x 10 ^ and radius p = 2.0 x 10 ®. 
By using the Matlab function eig to compute all eigenvalues of the test problem in 
dense format, we find that there are 141 eigenvalues within P. 

Define 


max_r = max r^, (4.1) 

1<1<S 

where Vi are the residual norms given by || Ax^ — Aii?Xi||/(|| Ax^||-|- ||i?Xi||) and (Ai,Xi) 
are the filtered eigenpairs. In our CIQZ method, we stop the iteration process when: 
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Fig. 4.1. The convergence behavior of CIQZ for test problem BFW782. 


(i) the numbers of filtered eigenpairs in two consecutive iterations are the same, and 

(ii) max_r in current iteration is less than the prescribed tolerance e (step 21 in Algo¬ 
rithm 3). 

The left picture in Fig 14.11 depicts the numbers of filtered eigenpairs in ten iter¬ 
ations. Recall that the filtered eigenpairs are the ones whose residual norms are less 
than the test tolerance rj. In the experiments, we take rj = 1.0 x 10“^. We see that in 
the first iteration there are 28 approximate eigenpairs whose residual norms are less 
than rj. But from the second iteration to the last, the number of filtered eigenvalues 
in each iteration is 141, which is exactly the number of eigenvalues inside F. 

The right picture in Fig 14.11 shows the maximum of the residual norms of filtered 
eigenvalues, i.e., max_r (cf. ()4.1ll l. in each iteration. From the left picture, we know 
that the number of filtered eigenvalues attains the one of eigenvalues inside F starts 
from the second iteration. Therefore, we plot max_r starting from the second iteration 
to the 10th. We see that maxjr decreases monotonically and dramatically from the 
second iteration to the fourth, maintains at almost the same level in the next three 
iterations, and rebounds from the eighth iteration. 

Experiment 4.2 This experiment is devoted to showing the numerical performance 
of our CIQZ. We compare the CIQZ method with Matlab built-in function eig 
and the block CIRR method (Block_CIRR). In [23], the authors addressed some 
implementation problems of the block CIRR method, including the selection of the 
size of the starting vectors and iterative refinement scheme. In the experiment, as 
for the block CIRR method, we use the version proposed in [23]. Note that the 
dominate computation cost of both CIQZ and Block_CIRR comes from solving 
g = 16 generalized linear shifted systems of the form (13.171) . On the other hand, when 
using eig to compute the eigenvalues inside the target region, we have to first compute 
all eigenvalues in dense format and then select the target eigenvalues according to their 
coordinates. However, the matrices listed in Table l4T] are sparse. Therefore, for the 
sake of fairness, we compare the three methods only in terms of accuracy, and will 
not show the amount of CPU time taken by each method. 
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Table 4.2 

Comparison of eig, Block_CIRR, and CIQZ. 


No. 

7 

P 

s 

eig 

Block CIRR 

GIQZ 

1 

-5.0 X 10^ 

1.0 X 10® 

58 

2.57 X 10-^® 

3.60 X 10-° 

4.76 X 10-1® 

2 

-6.0 X 10^ 

2.0 X 10® 

141 

5.59 X 10-12 

1.45 X 10-® 

2.09 X 10-12 

3 

-5.0 X 10® 

3.0 X 10® 

143 

6.81 X 10-1° 

1.41 X 10-® 

5.85 X 10-1° 

4 

-1.0 X 10^ 

8.0 

72 

1.15 X 10-® 

5.23 X 10-® 

7.00 X 10-1° 

5 

-5.0 X 10^ 

3.0 X 10^ 

137 

4.94 X 10-^ 

9.77 X 10-® 

1.67 X 10-° 

6 

-5.0 

3.0 

208 

1.99 X 10-® 

— 

1.79 X 10-® 


Block_CIRR and CIQZ are contour-integral based eigensolvers, the common 
parameters /iq and g we take to be 20 and 5, respectively. In [23], Block_CIRR 
performs two iterative refinements, i.e., three iterations in total. For comparison, in 
the experiment, we set the convergence tolerance e = 1.0 x 10“^^ and maxJter = 3 
for our CIQZ method. As a result, the results computed by CIQZ and Block_CIRR 
will actually be those computed in the third iteration. 

We use max_r (cf. (|4.1ll l to measure the accuracy achieved by each of three 
methods. In Table HaI 7 and p represent the center and the radius of target circle 
r respectively, and s is the number of eigenvalues inside F. In the last three columns 
of Table 14.11 we display the maxjrs computed by all three methods for each of the 
six test problems. 

From Table 14.11 we see that for the two contour-integral based eigensolvers, 
Block_CIRR and CIQZ, the latter outperforms the former in all six test problems. 
Especially, as for the problem 7, whose matrices are ill-conditioned, Block_CIRR 
fails to compute the desired eigenpairs. Therefore, our CIQZ method is more accu¬ 
rate and reliable than Block_CIRR. When it comes to the comparison of Matlab 
function eig and CIQZ, it is shown that the results computed by the two methods 
agree almost the same digits of accuracy for the first three test problems; while our 
CIQZ method is more accurate than eig by around two digits of accuracy in the last 
three problems, whose matrices are ill-conditioned. We should point out that in the 
experiment our CIQZ method just runs three iterations, it may obtain more accurate 
results if it performs more iterations. 

5. Conclusions. In this paper, we present a new contour-integral based method 
for computing the eigenpairs inside a given region. Our method is based on the CIRR 
method. The main difference between the original CIRR method and our CIQZ 
method is the way to extract the desired eigenpairs. We establish the mathematical 
framework and address some implementation issues for our new method. Numerical 
experiments show that our method is reliable and accurate. 
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